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We present a method that allows to distinguish between nearly periodic and strictly periodic time series. To 
this purpose, we employ a conservative criterion for periodicity, namely that the time series can be interpolated 
by a periodic function whose local extrema are also present in the time series. Our method is intended for the 
analysis of time series generated by deterministic time-continuous dynamical systems, where it can help telling 
periodic dynamics from chaotic or transient ones. We empirically investigate our method’s performance and 
compare it to an approach based on marker events (or Poincare sections). We demonstrate that our method 
is capable of detecting small deviations from periodicity and outperforms the marker-event-based approach in 
typical situations. Our method requires no adjustment of parameters to the individual time series, yields the 
period length with a precision that exceeds the sampling rate, and its runtime grows asymptotically linear 
with the length of the time series. 


Classifying the dynamics of a system plays an 
important role in its understanding. One of the 
main classes of dynamics are periodic ones and 
as a consequence deciding about periodicity is 
important in various scientific fields. While sev¬ 
eral methods for this exist, there is a shortage of 
such that can discriminate between nearly peri¬ 
odic and strictly periodic dynamics, as such de¬ 
tails are usually lost under experimental condi¬ 
tions (at which these methods are aimed). How¬ 
ever, the latter does not hold for simulated deter¬ 
ministic systems, which are often used as mod¬ 
els to improve our understanding of real systems. 
With such systems in mind, we propose an effi¬ 
cient periodicity test for time series that is ca¬ 
pable of detecting small deviations from period¬ 
icity. We demonstrate our method’s capabilities 
and show that it outperforms existing methods 
in typical situations. While our approach is not 
aimed at analyzing experimental time series, as¬ 
pects of it may enhance existing or inspire new 
techniques for this purpose. 


I. INTRODUCTION 

Two central distinctions in the classification of dynam¬ 
ical systems are that between a chaotic dynamics and a 
regular one as well as that between a transient dynamics 
and a stable on^^H^. As all periodic dynamics are neither 
chaotic nor transient, being able to decide that a dynam¬ 
ics is periodic is hence a valuable asset when analyzing 
dynamical systems. We here consider this issue for the 
case of simulated time-continuous deterministic systems. 

A plethora of methods to identify periodicities in 
time series have been proposed in the past, which can 
roughly be divided into four approaches: Periodogram- or 
autocorrelation-based methods^lhAJ search and evaluate 
local maxima in the frequency domain or in the auto¬ 
correlation function. Epoch-folding techniqueJ^^HIll ^j-e 


based on sorting the time series’ values into bins ac¬ 
cording to their phase for a presumed period length 
and finding the optimal period length according to some 
test statistics. Recently, approaches have been proposed 
that search for repeating patterns, i.e., subsequences of 
symbolJi^®^. Finally, Ref. [111 proposed to search for lo¬ 
cal minima of the mean fluctuation function. 

With exception for the latter, neither of these meth¬ 
ods aims at detecting small deviations from periodicity, 
because under experimental conditions, for which these 
methods were designed, such details can rarely be cap¬ 
tured due to noise and other confounding factors. For 
example, there may be no indicative difference between 
the Fourier spectra of chaotic and similar periodic sig- 
nal^i^. In contrast, time series from simulated systems 
may reflect small deviations from periodicity, which can 
be of considerable interest in their analysis. Moreover, 
unless the period length is a multiple of the sampling 
interval, all of the above methods either do not provide 
a straight forwa rd way to decide whether the time series 
is periodicliiElHIll^ e.g., via a critical value of some test 
statistics, or test again st the null hypothesis that the time 
series is randorrP l^° l ^^l ^, and thus do not distinguish be¬ 
tween regular and chaotic or stable and transient time 
series. 

Using embedding technique^^nUi] or time series of an 
observable and its temporal derivativ^^HM]^ available, 
one can reconstruct the dynamics’s phase space. Its struc¬ 
ture and other properti es may then form the basis of 
several kinds of analyse^^^^lHllI _ some of which require 
a careful selection of parameters and pose a lot of pit- 
falls -, but only in case of sinusoidal time series is there 
a criterion for periodicitj^^. 

If the dynamics of a system is known, methods that 
go beyond the analysis of observables become available: 
The theory of dynamical systems allows for an analyt¬ 
ical classification, e.g., by analyzing the stable and un¬ 
stable manifolds o f fix points and separatrices between 
basins of attractiorP^. Another approach (often used for 
continuation problems) is finding a piecewise polynomial 
approximation of a periodic orbit by numerically solv- 
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ing the respective nonlinear equation systerrPSl. These 
approaches may, however, become unfeasible for com¬ 
plicated and high-dimensional dynamics. A method that 
does not suffer from these problems is computing the sign 
of the largest Lyapunov exponent!^. Here the difficulty 
arises that a zero sign (signifying a regular dynamics) 
can only form the null hypothesis and can thus only be 
rejected and not accepted. 

We here propose a method that allows to decide 
whether a time series complies with a comparably con¬ 
servative periodicity criterion and that is thus capable 
of detecting small deviations from periodicity, such as a 
slowly rising amplitude or period length. We first moti¬ 
vate and define this criterion (Sec. IIA) and then describe 
our approach to tell whether it is fulfilled for a given 
period length (Sec. IIB), for any period length from a 


IICI and finally for any reasonable 


IIDI. We then extend our method 


small interval (Sec. 
period length (Sec. 
to account for small, bounded errors, such as nu merica l 
ones (Sec. II E[ ). Using archet ypal n oise-free (Sec. |III A ) 
and noise-contaminated (Sec. Ill B) test cases, we inves¬ 
tigate our method’s performance and compare it to an 
alternative approach based on marker events. Finally, we 
demonstrate that our method can help classify ing th e dy¬ 
namics of simulated dynamical systems (Sec. Ill CI and 
draw our conclusions (Sec. IVI. 


II. METHOD 

Some general remarks on notation: Lowercase italic let¬ 
ters denote integers; fraktur denotes fractions; Greek let¬ 
ters denote real numbers. Uppercase letters denote re¬ 
spectively valued functions. Blackboard-bold letters de¬ 
note sets. The binary modulo operator ranges between 
addition and multiplication in the order of operations. 

The source code of an implementation of our method 
in C as well as of a standalone program and a Python 
module building upon it is freely availabl^^. 


A. General approach and criterion 

Before delving into the mathematical details, we give 
a brief overview and illustration over our approach: We 
consider a non-constant time series 0 : {0,..., n — 1} —>■ 
M of length n > 3 , assuming without loss of generality 
a sampling rate of I. We want to decide about the ex¬ 
istence of an extension <I> of 0 to the interval [0, n — 1 ] 
that is periodic with a period length r (which we will 
initially consider to be fixed) and all of whose local ex¬ 
trema are captured by 0 . The latter criterion serves to 
avoid strongly oscillating solutions for <I>, which could al¬ 
most always be found. The function $ can be understood 
as the signal underlying 0 (see Fig. [^. 

Following the ansatz of epoch-folding techniqueJi^, 
we fold all periods into one by mapping (t, 0(t)) to 


(t mod r, 0(t)), thus obtaining an approximation of a sin¬ 
gle presumed period of $, which we refer to as foldation, 
T-foldation, or 0 .^ in the following (see Fig. [^. We then 
check whether the foldation’s number of local extrema 
complies with that of 0 . As an example, we show in Fig.[^ 
T-foldations from four time series, all of which have ca. 
two local extrema per r time units. We observe that small 
deviations of the period length or from periodicity suf¬ 
fice to drastically increase the foldation’s number of local 
extrema. 

As the exact definition of an individual local extremum 
is not relevant for our method, we only define the num¬ 
ber of local extrema of some function S. To this pur¬ 
pose we employ the longest zigzagging sequence with 
points from S’s graph ((70, S(7o)),..., (yfc+i, 5(7fc+i)), 
see Fig. [^: 

Definition 1 . We define the number of local ex¬ 
trema £’(S) of a function S : D —>■ K with D C M as 
the largest k such that there is an increasing sequence 
7 o,... , 7 fc-i-i G D for which 0 sgn(E;(7j) - E:(7 j_i)) 
sgn(S(7i+i) - S(7i)) 0 for all i e {1,..., /c}. 

Next, we define a correction term for extrema of some 
function S that 0 does not capture, but could capture if 
it began earlier or lasted longer, i.e., extrema that 0 fails 
to capture due to its finiteness and not due to S strongly 
oscillating (also see Fig.[^. 

Definition 2 . We say that S : D —>■ M with min(]D)) = 0 , 
max(D) = n — 1 , and { 0 , ...,n— 1 } C D has unde¬ 
tectable initial extrema if there is a A S D H [ 0 , 1 ] such 
that sgn(S(l) — S( 0 )) sgn(S(A) — S( 0 )) 7^ 0 . Analo¬ 
gously, we define undetectable final extrema. Finally, to 
correct for undetectable extrema, we define (^(S) to be 2 
if S has both undetectable initial and final extrema, 0 if 
it has neither, and I otherwise. 

Finally, we define the property, we are testing for: 

Definition 3 . We say that a time series 0 : 

{0,..., n — 1} —)■ M complies with a period length r iff an 
extension $ : [0, n — 1] —)■ K 0/ 0 exists such that $(t) = 
$(t mod t) yt e [0,n — 1] and U($) — (^(d)) = A(0). 

From this definition, it follows that 0 complies with every 
period length larger than n — 1 . Also, 0 complies with 
no period length shorter than 2, if n is sufficiently large. 


B. Deciding about periodicity with a given period length 

The following theorem (see Appendix|^for a proof) al¬ 
lows us to decide as to whether 0 complies with a period 
length r by counting the local extrema of an approxima¬ 
tion Qt of all of $ from a foldation (see Fig. Ilj). 

Theorem 1. For a given t, let D,. := 

{a G [ 0 , n — I] I G { 0 ,..., n — 1 } : i mod t = a mod r}. 
For a G Dt, let la G { 0 ,..., n — 1 } be defined such that 
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FIG. 1. Illustration of several concepts on the example of a time series 0 (green squares) defined on {0,..., 16} using a period 
length of r « 6.2 for folding. <I> (red line) is a r-periodic extension of 0. Except for an undetectable final extremum (see below), 
all of "F’s local extrema are captured by 0 and thus 0 complies with a period length of r and meets the criteria of Def. 
respectively. The arrow illustrates the process of folding ((t, 0(t)) —>■ {t mod r, 0(t))) for t = 9. The r-foldation 0r is marked 
by small dots; 0 is marked by circles. 7 o,. •. , 7 io indicate the values of one longest sequence 7 for counting the extrema of 0 
via Def. and 0 t have an undetectable final extremum at t « 15.7, which could be captured by 0 if it lasted until f = 17 
and contained, e.g., the point marked by the dashed square. <!> and 0 t have no undetectable initial extrema. 
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FIG. 2. Approximations of one presumed period of length r = 
27r (r-foldation) from exemplary time series 0 of length n = 
100. From left to right: 0(t) = sin(t); 0(t) = sin(1.003t); 
0(t) = (1 -I- 0.0011) sin(f); 0(t) = sin((l-|-0.000051) t). The 
numbers in the top right indicate the approximation’s number 
of local extrema. 


la mod T = a mod r. Finally define 0,- : D,- —>■ K 
via 0t(q:) := Q{la)- Then 0 complies with a pe¬ 
riod length t, if and only if 0 ,- is well-defined and 
E{er)-c{er) = E{e). 

0 r can only be ill-defined if there are € 

{ 0 ,..., n — 1} with i j and i mod t = j mod r, which 
in turn happens, iff r S ]F„, where F„ C Q is the set 
of fractions with a numerator smaller than n, i.e., the 
element-wise inverse of the n-th Farey sequence Fji (ig¬ 
noring the latter’s restriction to [ 0 , 1 ]). 

Instead of determining E{Qr), it often suffices to re¬ 
gard the number of extrema of the foldation 0 ,-, - i.e., the 
restriction of 0 t- to [ 0 , r) nD,- - and use it to extrapolate 
a lower bound for E{Qr), namely: 


-( 0 ) :=E{Qr)['^\+B{er){['^\-l)<E{e,), ( 1 ) 


period length r, Lt{Q) is often much larger than Eifd) 
(see also Fig. [^. In this case it suffices to regard a few 
values of the foldation to reject that 0 complies with a 
period length r. 

To determine the values of 0 t, functions with the fol¬ 
lowing property are useful: 

Definition 4. We say that a function I sorts the first b 
integers modulo r iff it maps { 0 ,..., 6 — 1 } to a permu¬ 
tation of itself and 

I{i) mod T < I{i 1) mod r Vf S {0,... , & — 2}. 

The function / that sorts the first n integers modulo r 
is unique and thus strictly monotonically increasing, if 
and only if r ^ F„. In this case, the values of 0 o / are 
identical to the values of 0 and assumed in the same or¬ 
der and, in particular, £1(0 o I) = E[Qr)- Moreover, we 
only need {0(/(O)),..., Q{I{n — 1))}, n, and I~^{n — 1) 
to determine the value sequence of Qr and thus E{Qr)- 
With additional knowledge of £“^(1) and I~^{n — 2) we 
can also determine C{Qr)- As knowing £l(0r) and C(0 t) 
(as well as £ 1 ( 0 )) suffices to decide whether 0 complies 
with a period length r (Th. [^, the latter only depends 
on I, with no further explicit dependence on r. 


C. Deciding about periodicity with period lengths within 
small intervals 

If we know r, we can easily find a function that sorts 
the first n integers modulo t by sorting. For ranges of r, 
we require the following theorem (see Appendix for a 
proof): 


where B{Qt) € { 0 , 1 , 2 } is the number of extrema that 
are so close to multiples of r that they are not accounted 
for by £l(0r). Moreover, if 0 does not comply with a 


Theorem 2. Let 1 < t < n and let | be the largest 

and j be the smallest reduced fraction from F„ such that 
I < T < j. Define £„,r(f) ■= ir mod {p-\-r) for i G 
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{0, ...,p + r}. Then In^r sorts the first p + r integers 
modulo T. Moreover, /„ t increases strictly monotonically 
on {0,... ,n - 1} ifT< 

It follows that, for two successive elements a and b of 
F„ with a < b, is the same for all r S (a, b]. There¬ 
fore, if 0 complies with one period length r in (a, b), it 
complies with all period lengths in that interval. More¬ 
over, if 0 complies with a period length b, it also com¬ 
plies with all period lengths in (a, b) and thus it suffices 
to investigate open intervals of adjacent elements of F„ to 
decide about the periodicity of 0. As a consequence, all 
we need to do to check whether 0 complies with period 
lengths in (o, b) is to count the extrema of one foldation. 

In another consequence, Th. allows us to iterate over 
the values of In,T and thus of 0 o without explicitly 
calculating In,T{i) for all* G {0,..., n — 1}. As we men¬ 
tioned in the previous subsection, a few such iterations 
often suffice to reject that 0 complies with a period r. 


D. Deciding about periodicity with an unknown period 
length 


If we have little prior knowledge or constraints on the 
value of a possible period length t, testing every possi¬ 
ble interval of adjacent elements of F„ will usually be 
unfeasible because |F„| = O(n^) as n —>■ ocP. To avoid 
this, we can make use of the following: First, if 0 com¬ 
plies with a period length r, so does any segment of 0. 
Second, if the length n of such a segment is sufficiently 
smaller than n, the interval of adjacent elements of F^ 
in which a given r lies is larger than the analogous in¬ 
terval for F„. Third, for each two adjacent elements on 
some level of the Stern-Brocot tree^S, there is some o 
such that they are adjacent elements of Fq- Combining 
these three facts, we can make a nested-interval search 
on the Stern-Brocot tree for a r with which 0 complies, 
as follows: 

Method: To test whether 0 complies with some pe¬ 
riod length in 

ments of F^ for some m < n: 


with I and ^ being successive ele- 


(1) If p-|-r > n: Check whether 0 complies with period 

lengths from (using Th. 1). 

If yes, the test is positive. If no, the test is negative. 

(2) If p -I- r < n, check whether some segment of 0 of 
length n := p + r complies with period lengths in 

(|, j) (using Th. j^. If j = g, this check can be 

skipped as it is always positive. 

If no, the test is negative. If yes, continue. 


( 3 ) 


Use this method to test whether 0 complies with 
some period length in ( ^. 


If yes, the test is positive. If no, continue. 







7 18 n 

2 5 3 


FIG. 3. Illustration of a possible realization of the nested- 
interval search for a time series with 18 < n < 25 that com¬ 
plies with period lengths in (^, . Each level corresponds 

to one level of the Stern-Brocot tree. The numbers below the 
intervals indicate the length of the segment checked in step 
or[^of the procedure for the respective interval, i.e., n or n, 
respectively. “A” indicates an automatic pass due to ^ be¬ 
ing the interval’s right border. Checkmarks indicate whether 
this check was positive or negative; only displayed intervals 
were checked at all. Green, straight lines indicate that the 
test found 0 to comply with some period length in the re¬ 
spective interval; red, zigzagged lines indicate that it did not. 
The check on the interval (f, |] is an example for a false 
positive result in step[^ Intervals are not displayed to scale. 


(4) Use this method to test whether 0 complies with 
some period length in p • 

If yes, the test is positive. If no, the test is negative. 

Without any prior knowledge or constraints on the value 
of a possible period length r, the above test can be ap¬ 
plied to (j, ^], i.e., with p = 2, q = r = 1, and s = 0, 
neglecting all intervals whose smallest value exceeds some 
Tinax < n—1. A realization of this procedure is illustrated 
in Fig. 

Some remarks on the implementation: 

• If the check in step yielded a false positive re¬ 
sult (such as for the interval (j, |] in Fig. this 
would not affect the total outcome of the tesrj as all 
subintervals will be tested again (in step or at 
a higher recursion level. (Some segment of 0 may 
comply with a period length r while 0 itself does 
not, anyway.) This can be employed to make an im¬ 
plementation more effective by using the number of 
extrema of the foldation of the respective segment 
to extrapolate the number of extrema of 0r and 
comparing it to U(0) (see Eq. [^. This way, cal¬ 
culating or storing the number of local extrema of 
each segment can be avoided. 

• We empirically found that the runtime was in¬ 
creased if we used a random segment of length p+r 
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instead of the first one in step 


• As the interval for some m < n re¬ 

sides at the (to — l)-st level of the employed branch 
of the Stern-Brocot tree, a fully recursive imple¬ 
mentation can result in a problematically high re¬ 
cursion depth for large n and r„iax- This can be 
avoided by applying the above test to the inter¬ 


vals 


(2 3] (3 4] f [Tmaxl-l 


C'^maxl 

i 


in a non- 


recursive manner. 


• The above method can also be used to hnd the 
shortest period length that 0 complies with. For 
obtaining an estimate with an error margin for the 
shortest period length of the underlying process, 
one usually would want to know the first maximal 
interval such that 0 complies with every period 
length in that interval. To avoid underestimating 
the margin, one needs to consider that if c is the 
successor of the successor of a in ]F„, 0 may com¬ 
ply with all period lengths in (a, c]. Therefore, after 
finding an interval such that 0 complies with ev¬ 
ery period length in that interval, the next interval 
of adjacent elements of F„ needs to be checked as 
well. 

In the following, we consider and use an implementa¬ 
tion that takes all of the above into account. The asymp¬ 
totic runtime of our method scales linear with n, if 0 
complies with a period length t <^n and quadratic with 
Tmax otherwise (see Appendix [C|). 


E. Accounting for small errors 


In this subsection, we describe a simple expansion of 
our method that is capable of allowing for errors of the 
time series 0 that are bounded and small in comparison 
to 0’s features of interest. Such errors could originate 
from small numerical inaccuracies, e.g., of a solver for 
differential equations, or small measurement errors. We 
handle such errors by treating two values of 0 that are 
less apart than a given error allowance a as identical 
for the purpose of determining local extrema. More pre¬ 
cisely, we weaken our criterion (Def. by changing Defs. 
[Hand[2] as follows: 

• In Def. 12 we additionally require that 
|S( 7 i+i) - S( 7 i)| > cr for all ie {0,...,k}. 


• In Def. if i and j with i < j are the smallest 
elements of {0,..., n — 1} such that |5(j) — S(*)| > 
CT, we say that S has undetectable initial extrema if 
there exist A,/z S D n [0,j] such that \fi — X\ > a 
and sgn(^ - A) ^ sgn(S(j) - E{i)). 


For cr = 0, the adapted Def.j^is equivalent to the original 
one, and so is the adapted Def. unless 0(0) = 0(1) or 


to 


0(n — 2) — 0(n — 1). The remainder of Sections IIA 
jllDj as well as the appendices analogously apply to the 
expanded method. 


207r 10® 10'‘ 10® 10® n.p, 
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10 "®® 

10® 10® lO'® 10® 10® 10® 10® 10"® 10® 10® 

n n 

FIG. 4. Smallest period length r < n — 1 with which the 
aperiodic time series (left) or (right) complies for 
f = 207r. Black points (n.p.) indicate that no such period 
length exists. 



III. APPLICATION 


A. Data with known periodicities or aperiodicities 


To investigate the performance of our test, respectively, 
we employed the following families of time series that 
deviate from periodic ones to an adjustable extent: 

Tenf : {0,..., n — 1} —1 

Tenf(t) := (l + e^) 

^enf ■ {0, . . . , n — 1} —>■ 

/ 27rt 

■^,nf{t) := cosf 

' T 

with e, f S M“''. The time series T^nf features a rising 
amplitude; the time series '^enf features a rising period 
length. The parameter e determines how quickly the am¬ 
plitude or period length, respectively, is rising, i.e., how 
strongly the time series deviates from a periodic one (at 
e = 0), whose period length is f. 

In Fig. we show the smallest period lengths r with 
which these time series comply depending on n and e for 
f = 207r. For both families of time series, we find that 
T « f for small n and e, more precisely for en® ^ 10®, 
i.e., for a relative amplitude/frequency change smaller 
than For higher n and e, the period length t ten- 
dentially increases with e and n, with r obtaining higher 
values in general for than for T^nf. Most time series 
with high n and e (roughly: en > 10"®, relative ampli¬ 
tude/frequency change bigger than 10"^) do not comply 
with any period length smaller than n — 1 and if they 
do, T is close to n. For f = 27r, r exhibits comparable 
dependencies on n and e, however, with some anomalies 
due to aliasing effects (not shown). 

These results show that our test is capable of detecting 
small deviations from periodicity. However, there is a con- 
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FIG. 5. Left column: Smallest cf such that the aperiodic time 
series (top) or (bottom) does not comply with a 

period length less than 1.5 f as per our foldation-based test. 
Central column: Smallest tm such that (top) or 

(bottom) are aperiodic as per the marker-event test. Right 
column: Quotient of and et. For T, e^i is larger than et 
everywhere, except if Cm = 0 (see top middle panel). For Vl', 
Em is larger than Cf in the top left half of the diagram (except 
if Cm = 0) and smaller in the bottom right half. The dashed 
black line in the bottom right panel marks the parameters 
investigated in Fig. 

siderable set of parameters for which periodicity with a 
long period length is detected. Moreover, even time series 
with strong deviations from periodicity sometimes com¬ 
ply with period lengths close to n. This can be explained 
by the fact that a time series 0 can easily “accidentally” 
comply with a period length close to its length n, e.g., for 
0 to comply with period lengths from (n — 2, n — 1), it 
suffices that 0(n — 1) lies between 0(0) and 0(1). This 
demonstrates the importance of choosing r„iax properly. 

We now compare our test’s performance to a test based 
on Poincare sections or marker event^^, which we refer 
to as marker-event test. As marker events we employ on 
the one hand the upward zero crossings of a piecewise lin¬ 
ear interpolation of the time series and on the other hand 
the time series’s local maxima. We consider a time series 
periodic according to the marker-event test, if neither 
the distances of subsequent zero crossings nor the ampli¬ 
tudes of local maxima are significantly correlated with 
time (with a significance level of 0.05 as per Kendall’s 
correlation coefficient). To determine the period length 
with this test, we use the mean of the distances of sub¬ 
sequent zero crossings, and as its confidence interval, we 
use twice the standard error. 

As a first benchmark, we use the lowest deviation from 
periodicity e for which the test detects T^nf or 'i’enr, re¬ 
spectively to be aperiodic (Fig. [^. We find both tests to 
be more specific for higher n and f in general. However, 
the marker-event test detects even the purely sinusoidal 
Tonf = 4'onf to be aperiodic in some instances, which is 
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FIG. 6. Left: Relative size St of the confidence interval for the 
period length for our test for sinusoidal time series with period 
length r and n sampling points. (5f := with (o, b) being 
the first maximal interval such that the time series complies 
with every period length in that interval.) Genter: relative 
size Sra of the conhdence interval for the marker-event test, 
(dm := with ? being the standard error of the distances 
between adjacent upward zero crossings.) Right: Quotient be¬ 
tween the two. dm is larger than df in the top left half and 
smaller in the bottom right half. 


expected given the possibility of type I errors by the sta¬ 
tistical test. Apart from these cases, our test is generally 
capable of detecting smaller deviations from periodicity 
in the form of increasing amplitude (T) than the marker- 
event test (blue points in Fig. top right). For changes 
of the period length (\k), the marker-event test performs 
better for high r and small n, namely for n ^ 

As a second benchmark, we regarded the error margin 
of the estimate of the period length for sinusoidal time 
series (Fig. [^. For both tests, we find the margin to de¬ 
crease with T and n, however, the latter decrease is small 
for the marker-event test. For high t and low n, the er¬ 
ror margin is higher for our test (blue points in Fig. 
right), while it is higher for the marker-event test for 
low r and high n, more precisely for n ^ For 5305 
of the 90000 time series analyzed for Fig.[^ the actual pe¬ 
riod length (r) did not lie inside the marker-event test’s 
error margin, while it did always lie within the margin 
for our test. 

Our results show that our test outperforms the marker- 
event test for coarse sampling and a high number of data 
points as well as for rising amplitndes. Moreover, it has 
no false positives and can thus be regarded to be more 
robust. Note that the marker-event test used here was 
tailored to the investigated time series (by assuming one 
upward zero crossing and one local maximum per pe¬ 
riod) and to the types of deviations from periodicity (by 
assuming a rising amplitude or period length), while our 
test requires no comparable adjustment. Finally, if Tmax 
is not chosen too high, our test’s asymptotic run-time be¬ 
havior is better than that of the marker-event test, which 
is 0{nlog{n)) due to the correlation coefficienfP^. 
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FIG. 7. Bottom: Smallest error allowance cTmin for which sinu¬ 
soidal time series of length n with a period length r = 2007r 
that were contaminated with white noise from W([ 0 , v\) com¬ 
ply with a period length of r. The black solid line indicates 
Cmin = v- Arrows indicate the third-smallest distance be¬ 
tween adjacent values of the r-foldation. Top: Relative size 
of the confidence interval for the period length (cf. Fig. [m 
for n = 10® and an error allowance a — v. For n = Im, 
a comparable behavior was observed (with being gener¬ 
ally larger); for n — 10“^ and n = 10®, was constant for 
10~®® < V < 10“^. In both plots, the and percentile 
over 200 realizations of the noise are shown. 


B. Noisy data 

To investigate the impact of erroneous data on our 
test, we first apply it to sinusoidal time series that are 
contaminated with white noise from 77([0, iz]), with lA de¬ 
noting the uniform distribution. In the bottom part of 
Fig. [tI w e show the minimum error allowance CTmin (see 
Sec. llIE[ ) that needs to be made for such a time series 
to comply with the correct period length. We find that 
the noise does not affect the result up to a certain noise 
amplitude zzcrit- For up to roughly 61 /critj the noise’s im¬ 
pact is strongly realization-dependent. For higher noise 
amplitudes, CTmin is slightly smaller than the noise ampli¬ 
tude V. We explain these regimes as follows: For v ^ zzcrit, 
the distances between consecutive values of the foldation 
are larger than v and thus the noise cannot introduce ad¬ 
ditional local extrema. This is confirmed by the observa¬ 
tion that r'crit roughly corresponds to the third-smallest 
such distance (the smallest and second-smallest occur at 
local extrema, where changing the order of values does 
not affect the test’s outcome; see arrows in Fig. [^. For 
v ^ biZcrit, the probability that the noise did not intro¬ 
duce any additional local extrema to the foldation be¬ 
comes negligible and thus an error allowance of roughly 
V is needed for compensation. Moreover, we find that for 
small noise levels the accuracy of the identified period 
length is not affected (see the top part of Fig. [^. We 
conclude that ct = zz is an appropriate choice, given a 



FIG. 8 . Minimum noise level z/min such that ^4 ^ 2 . 4 , f 

(see text) contaminated with a white noise from W([—f, |]) 
is (wrongly) detected to be periodic with a period length of 
about T according to our test (with an error allowance a = v) 
or the marker-event test, respectively. 


known noise amplitude v. 

To evaluate our test’s robustness against noise and to 
compare it with the marker-event test, we employ '^enf 
contaminated with white noise ^([—|, |])- To exclude 
other factors that may influence the relative performance 
of the tests, we set n = for which « Cf for the 
uncontaminated time series (see the dashed line in the 
bottom right of Fig. [^. For these cases we found that 
Cm ~ Cf ~ ^('?) = 4.1 • with the coefficients of 

the latter being obtained by a fit. We chose e = Cfl(r) 
with ^ = 10. In Fig. we show the minimum noise level 
i^min for which our test (with an error allowance a = v) 
or the marker-event test fail to detect the aperiodicity 
of for the above conditions. For our method, we 

observe zzmin to slightly increase with f for low f , being 
higher than the amplitude of the uncontaminated time 
series. The time series that were detected to be periodic 
were constant with respect to the error allowance ct, i.e., 
max(0) — min(0) < ct. Around f « 27r, Vmin quickly de¬ 
creases to about 0.5, after which it decreases more slowly 
with some power law. We explain these two regimes as 
follows: If v (and thus ct) is smaller than some value be¬ 
tween 0.5 and 3, periodicity is detected when the error 
allowance prevents the test from detecting the uncon¬ 
taminated time series’s deviations from periodicity. If, 
however, v is larger, the noise dominates the original sig¬ 
nal and the original signal acts like a contamination (and 
thus the condition that the errors are small in compar¬ 
ison to the time series’ features of interest is not met 
anymore). Under an error allowance ct = zz, the noise 
is constant in the terms of our test and thus periodic. 
Thus, periodicity is only detected, when the noise am¬ 
plitude becomes so high that the uncontaminated time 
series’s influence on the noise becomes negligible. For the 
marker-event test, i^rnin mostly follows a power law, the 
main exception being a few cases, in which zzmm = 0. We 
made comparable observations for ^ = 100 and ^ = 1000. 
In general is higher for our test, which indicates that 
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FIG. 9. Minimum error allowance cTmin such that a time se¬ 
ries of length 10® generated by a model system of two coupled 
FitzHugh-Nagumo oscillators (see text) integrated with dif¬ 
ferent adaptive Runge-Kutta methods with a maximum rel¬ 
ative error p and a maximum step size of 1 complied with 
some period length smaller than 500. The solid line indicates 
Umin = p. Only cases with Umin yf 0 are depicted. The num¬ 
ber of cases with (Tmin = 0 were: Kutta: 37; classical: 153; 
Fehlberg: 112; Cash-Karp: 154; Dormand-Prince: 160 (all). 


it is less affected by noise, provided the error allowance a 
can be chosen to match the noise amplitude. 


C. Dynamical Systems 

To evaluate our method’s performance on the analysis 
of dynamical systems, we apply it to a deterministic sys¬ 
tem of two diffusively coupled FitzHugh-Nagumo oscil- 

latorpHm 

We employ a parameter range, in which this 
system exhibits several regimes of (periodic) mixed-mode 
oscillations (MMOs)P^ separated by chaotic windows. To 
describe these MMOs, we use the following notation: 
h^ih ^2 ■ ■ -I which indicates that one period consists of hi 
high-amplitude oscillations, followed by li low-amplitude 
oscillations, followed by ^2 high-amplitude oscillations, 
and so on. 

We chose an initial condition near the attractor, dis¬ 
carded transients and integrated this system’s dynamics 
for n = 10® time units, sampling each time unit, with 
several adaptive Runge-Kutta methods, namely Kutta’s 
3'^‘^-order method (using the midpoint method for error 
estimation)!^, the classical Runge-Kutta method (with 
step doubling used for error estimation), Fehlberg’s d**'- 
order method, the Cash-Karp method, and Dormand’s 
and Prince’s 8*^-order method - all as implemented in 
the GNU Scientihc Librarj!^. We applied our test to the 
temporal evolution of the hrst oscillator’s first dynamical 
variable (xi in Ref. I35[) . whose maximum absolute value 
was ca. 0.9. Due to the latter, the highest expected abso¬ 
lute integration error roughly corresponds to the relative 


integration accuracy p. 

For a first analysis, we chose a coupling strength k of 
0.17, for which the system exhibits a 1^ MMO, i.e., a pe¬ 
riodic dynamics (see inset of Fig. 10). We considered the 
test to be successful if the respective time series complied 
with some period length smaller than 500 - with the pe¬ 
riod length of the system’s dynamics being roughly 287. 
In Fig. we show the minimum error allowance CTmin 
needed to be made for the test to be successful depend¬ 
ing on the integration accuracy p. For Kutta’s 3''‘^-order 
method, (Tmin is mostly one order of magnitude smaller 
than p for 10“^ ^ p ^ 10“^ and only larger than p in 
one case. For all other integration methods, Umin is never 
larger than p and deviates from 0 for at most a few val¬ 
ues of p, which exhibit no discernible pattern (except for 
the Fehlberg method and p > 10“^). While a detailed 
investigation of this phenomenon and why it does not af¬ 
fect Kutta’s 3"'‘^-order method is beyond the scope of this 
study, we hypothesize that it can be explained as follows: 
The cases with (Tmin = 0 are due to the fact that the in¬ 
tegration error is not stochastic but systematic in nature 
and therefore likely to affect adjacent values of the folda- 
tion in a comparable way, thus not affecting their order. 
We hypothesize that this effect is diminished if the step 
size is frequently adapted, which leads to the seemingly 
random deviations of (Tmin from 0. From the above results 
and our results from the previous subsection, we conclude 
that the maximum expected absolute integration error is 
a good choice for a. 

We conclude this section with two applications of our 
test to identify possible periodicities in time series gener¬ 
ated by simulated dynamical systems. First, we investi¬ 
gate the coupling regimes of the aforementioned system 
of two diffusively coupled two-dimensional FitzHugh- 
Nagumo oscillators. To this purpose, we employ the 
shortest period length t as found by our test and the 
temporal distances p between subsequent high-amplitude 
oscillations - a marker-event-based observable. In Fig.fTOl 
we show the dependence of r and p on the coupling 
strength k. While both allow to separate the regimes of 
the primary MMOs (1°, 1^, ...), only t clearly discrim¬ 
inates between the regimes of secondary MMOs (2^, 3^, 
...; iH®, IH^U, ...; I2l3^ 12 i2i3 ^ ) xhis demon¬ 

strates that our method may provide complementary in¬ 
formation when analyzing dynamical regimes - in addi¬ 
tion to telling chaotic dynamics from periodic ones. 


Finally, we apply our method to time series generated 
by a small-world network of lOOOO diffusively coupled 
FitzHugh-Nagumo oscillator^^ (the system will be dis¬ 
cussed in detail elsewher^^. We evolved these systems 
with Fehlberg’s 4*’*'-order method with a relative error 
p = 10“® and employed the average of the first dynami¬ 
cal variable {x in Refs. [M] and SOI) as an observable. As 
the first dynamical variables roughly range between —0.4 
and 0.9, we expect the maximal absolute integration error 
to be roughly p. This system is of interest here because 
it is capable of exhibiting long, nearly periodic episodes, 
which eventually turn out to be a transient behavior. 
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coupling strength k 


FIG. 10. Top: Observed lengths of intervals r] between subse¬ 
quent high-amplitude oscillations (as identified by threshold 
crossings of a piecewise linear interpolation of the observed 
time series) depending on the coupling strength ac in a sys¬ 
tem consisting of two coupled FitzHugh-Nagumo oscillators 
(see text) integrated with Fehlberg’s d^^'-order method and a 
relative error p = 10“® for n = 10® time units (data already 
shown in Ref. EH. Bottom: Lowest period length with which 
the respective time series comply for an error allowance of 
o — 10“®. The occasional single dots (e.g., between the 1^ and 
the 1^1^ regime) occur when a sampling point falls within one 
of the tiny periodic windows within the small chaotic windows 
that separate the larger regimes. Labels: Mixed-mode oscilla¬ 
tions corresponding to selected regimes. Inset: Excerpt of the 
time series for k = 0.17 (1^ MMO). 


In the top of Fig. El we show a time series contain¬ 
ing such an episode, whose aperiodicity becomes evident 
through its finiteness. For 1500 ^ t ^ 15000 this time 
series appears periodic with a period length of roughly 
650. Had our observation ended at t = 15000, visual in¬ 
spection might thus have led us to the false conclusion 
that the dynamics might have become a stable, periodic 
one. Exemplarily applying our test to this time series for 
the t G [3000,13000], i.e., to the center of the periodic 
episode, we find that this excerpt does indeed not com¬ 
ply with any period length t G [2,10000] for any error 
allowance cr smaller than 10“^'^. In the bottom of Fig. 11 


we show an excerpt of another episode, which continued 
for at least another 1.5 • 10® time units and which we 
thus hypothesize to be periodic. Applying our test to 
this time series for t G [3000,13000], we find that this 
excerpt does comply with period lengths around 388 for 
error allowances 10“® ^ cr ^ 10“^'^. We obtain the same 
result when applying the test to the time series for the 
following 1.5 TO® time units, which affirms our hypothesis 
that this episode is actually periodic. The marker-event 
test finds both (de-meaned) time series to be periodic 
and thus in particular fails to detect the aperiodicity of 
the first time series. In all cases, we obtain comparable 
results for other, similar intervals. 


IV. CONCLUSIONS 


We proposed a method to test whether, for a given time 
series, there is a periodic function that interpolates it and 
whose local extrema are captured by the time series. Due 
to the conservativeness of the criterion, our method is 
highly specific and capable of detecting even small de¬ 
viations from periodicity. Moreover, our approach yields 
an interval of possible period lengths that is usually nar¬ 
row in comparison to the sampling rate and allows for a 
precise reconstruction of one period of the observable of 
the time series. We found that, in typical situations, our 
method outperforms an alternative, marker-event-based 
test in terms of specificity, precision of the detected pe¬ 
riod length, and robustness - even though this test was 
tailored to the investigated time series. By applying it to 
two typical problems, we also demonstrated our method’s 
usefulness for the analysis of simulated time-continuous 
dynamical systems. 

The first parameter that needs to be chosen for our 
method is the maximal accepted period length. Choos¬ 
ing it too high may impair the specificity of the method 
and its runtime if the time series is not periodic with 
a short period length. The second parameter is the er¬ 
ror allowance, which can be straightforwardly chosen in 
the case of a simulated system via the error of the in¬ 
tegration method. The individual features of the time 
series do not affect the choice of either of these parame¬ 
ters - while they have to be taken into account for many 
other methods, e.g., when choosing the marker-events for 
a marker-event-based approach. 

For high-dimensional systems, the main computational 
challenge (in terms of both, runtime and stability) re¬ 
quired to apply our method is evolving the system’s dy¬ 
namics to generate a time series. Thus, if one has already 
performed the latter, our method can be applied with lit¬ 
tle effort—in contrast to commonly used techniques such 
as the maximum Lyapunov exponent or numerically find¬ 
ing an approximation of a periodic orbit. The latter is 
often used in continuation method^^, whom our method 
may assist by providing accurate starting values, namely 
a precise estimate of the period length as well as an ap¬ 
proximation of the orbit by folding all dynamical vari¬ 
ables. 

It is essential for our approach that the period length 
does not vary considerably and that there are no phase 
jumps or drifts - a requirement that is almost only ful¬ 
filled by deterministic and stationary systems. Moreover, 
our test requires errors to be small and bounded with 
respect to features of interest. These requirements are 
rarely met by real systems or experimental observations, 
respectively, and thus we cannot expect that our test in 
its present form will find application for experimental 
data. Nonetheless, parts of our approach, in particular 
our findings on the arithmetics of folding (Th. and 
the nested-interval approach, may enhance existing or 
inspire new epoch-folding techniques. Moreover, a varia¬ 
tion of our approach may be applicable to the analysis of 
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FIG. 11. Two exemplary time series of the mean of the first dynamical variable of 10000 FitzHugh-Nagumo oscillators that 
were diffusively coupled on a small-world topology. 



those systems that can be approximated as deterministic 
and stationary for a sufficiently long time with respect 
to the Mi'iod length. One example for such systems are 
pulsarP2, where the additional knowledge of the observ¬ 
able’s temporal derivative may be employed to make our 
approach applicable t o peri od lengths that are smaller 
than the sampling 
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Appendix A: Proof of Theorem 

Theorem 1. For a given t, let D,- := 

{a S [0,n — 1] I S {0 ,... ,n — 1} : i mod t = a mod r} 
For a € Dt, let la S {0,..., n — 1} be defined such that 
la mod T = a mod r. Finally define 0,- : —)■ K 

via Qt{cx) ■= 0(^c()- Then 0 complies with a pe¬ 
riod length t, if and only if 0,- is well-defined and 
E{er)-cier) = E{e). 

Proof. We first note that if some function r 2 extends 
some other function Fi, we have A(ri) — C'(ri) < 

A(r2)-C'(r2). 

If 0 complies with a period length r, there exists a d* 
that extends 0 to [0, n — 1] and with <l>(t) = mod r). 
We therefore have: 

0t-(q;) = Q{la) = ^{la) = mod r) 

= $(a mod r) = $(q;) Vo G Dt-, 

and thus <f> also is an extension of Or- Because of this, of 
Or extending 0, and of C(0) = 0, we have 

E(0) = E(0)-C(0) < E(0r)-C(0r) < E(<^)-C(^). 

Therefore, £!($)-C'(d>) = E(0) yields E(0r)-C(0r) = 
E(0). For 0r to be ill-defined, la must be ill-defined and 


thus, there need to be i,j G {0,... ,n — 1} with i mod 
r = j mod T and 0{i) fi 0{j). From this, it directly 
follows that no extension of 0 can be r-periodic. 

To show the other direction of the equivalence, we con¬ 
struct an extension $ of 07 - to [0, n — 1] by piecewise lin¬ 
ear interpolation, for which we thus have A($) — (7(4)) = 
E{0r) — C(0r). If E{0r) — C{0r) = E{0), we thus also 
have E{^) — C'($) = E{0). Furthermore, 

$(a) = 0ria) = 0{la) = 0ila mod r) 

= 0t(q: mod r) = $(a mod r) Va G Dt-, 

and due to piecewise linear interpolation, also $(t) = 
mod r) Vt G [0,n—I]. Thus, the requirements of 
Def. [3] are fulfilled. □ 


Appendix B: Proof of Theorem 

Theorem 2. Let 1 < t < n and let | be the largest 

and J be the smallest reduced fraction from F„ such that 
- < T < V Define ■= ir mod {p -\-r) for i G 

{0, ... ,p -\- r}. Then In^r sorts the first p -\- r integers 
modulo T. Moreover, E r increases strictly monotonically 
on {0,... ,n - 1} ifT< fi 

Proof. Since F„ is the element-wise inverse of the n-th 
Farey sequence F„, we directly get from the theory of 
Farey sequence^^ 


ps -\-1 = rq 


(Bl) 


Thus p and r are coprime and consequently so are r and 
p -\- r. From this we get that for eA\ i fi j] i,j<p-\- r: 

In.rfi) = ir mod {p -\- r) fi jr mod {p -{- r) = 

i.e., that is a bijection on {0, ... ,p -\-r -1^. 

From the definition of p, q, r, and s, we get: 


—p mod T = qr — p> 0 (B2) 

r mod T = r — ST > 0. (B3) 


Now, let fc G {0,..., p -b r — 1} and 


dk ■= 


kr — kr mod (p -\- r) 
p-\-r 


kr 

p-\-r 


<r-I. (B4) 
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As dk and k — dk are monotonically increasing with fc, we 
can make the following estimate for k — dk by inserting 
the lowest and highest value for k: 


0 < k — dk < {p + r — 1) — {r — 1) = p. (B5) 


Equations B2 to |B5| allow us to show the following in¬ 
equality: 


(fc — dk){r mod r) -I- dk{—p mod r) 

. . X , . X . X E] 

< p(r — ST) + r(qT — p) = [rq — ps)T = r. 

This also gives us that both, {k — dk){r mod t) and 
dk{—p mod r), are smaller than r. Using 

m S N; /3 G M; m(/3 mod r) < r 
mP mod r = to(/ 3 mod r), 

we can thus write: 

(fc — dk)r mod t = {k — dk){r mod r) (B7) 
dk ■ (—p) rnod t = dk{—p mod t). (B8 ) 

Finally, we can write: 

In,T{k) mod T 
= [kr mod [p + r)] mod t 

® [kr — dk{p + r)] mod T 
= [{k - dk)r + dk ■ (-p)] mod r 

= [{k — dk)r mod t -\- dk ■ (—p) mod r] mod r 

— dk){r mod r) -|- dk{—p mod r)] mod r 
^ (fc — dfe)(r mod r)-I-p mod r). 

Since k — dk and dk are both monotonically increasing 
with k and at least one of them increases if k is increased 
by 1, In,T{k) mod r is strictly monotonically increasing, 
unless T = - and thus r mod r = 0, in which case it is 
only weakly monotonically increasing. □ 


Appendix C: Asymptotic runtime behavior 

We first estimate the behavior of the average runtime, 
if 2 < r ^ n is the shortest period length that 0 com¬ 
plies with. To this purpose we employ the following facts, 
approximations and assumptions: 

(A) The runtime of the checks performed in steps 
and of the algorithm is approximately cn or ch, 
respectively, with some constant c. This is based 
on the assumption that run-time reductions due 
to aborting the counting of extrema early because 
their number already suffices to reject periodicity 
can be accounted for by a constant factor (which is 
already incorporated in c) . We approximate every¬ 
thing except these checks to have a runtime of 0. 


(B) We assume that, if 0 does not comply with any pe¬ 
riod length in the check in step 2 is positive 

with a probability c/> and that this is independent of 
other results. We further assume that (j) < ^, which 
we could confirm empirically for exemplary time se¬ 
ries without often repeating values, i.e., for which 
there were no pairwise different ti,... ,ty such that 
0 (ti) = ... = 0(U) and ^ > 0. 


(C) We assume that, at some level of the binary search, 
T is equally likely to be in the left or right branch 
of the Stern-Brocot tree and this is independent of 
other results. 


(D) On any level j of the Stern-Brocot tree, we approx¬ 
imate that, for each interval, fi corresponds to the 
average value of n over all intervals on this level. We 
denote this average value by fij. This approxima¬ 
tion is based on the assumption that it is essentially 
at random which intervals are investigated. 

(E) Let h be the level at which the smallest interval 
containing r involved in the search resides. Then 
we assume that no higher level than h is involved 
in the search. 


(F) One can neglect the additional checks that are 
made to ensure that a maximal interval is found 
such that 0 complies with all period lengths in that 
interval. 


Going by these assumptions, we now first estimate the 
number aj of checks performed on the j-th level (of the 
employed branch of the Stern-Brocot tree): At each level 
between 1 and h, we perform one check for the interval 
containing r. At each level between 2 and h an additional 
check is performed with a probability of ^ (if r is in the 
right branch; see assumption]^. Each of these checks has 
a probability of (j) to cause two additional checks an the 
next level, each of which in turn has a probability of (j) 
to cause two additional checks on the next level and so 
forth (assumption]^. We therefore obtain on average: 


J-2 


— 1 + 


^(2<^)*<1+- 


0 ifj = l 
I else 


oo .. .. 

y"(2(j)y = - + -=: d. 

2 ^ 1-26 
i=0 ^ 


i=l 


Now, let the numerators on some level j of the Stern- 
Brocot tree be zi,..., Zyj. Then the numerators on the 
next level are zi,zi + Z 2 , Z 2 ,..., Zyj-i + Zw, z^ and thus 
(by approximation ]d]) : 


W— 1 


- Ui + 2 E ■ 


z=2 


" w-1 

3 I 321 -b 6 E 


2 1 











12 


Using this, we obtain for the runtime gj of each check 
performed on the j-th level (by approximation [A| : 

gj = chh-i ( - j Vj e {1,..., - 1} and gh = cn. 


As fih-i < n , we can thus estimate gj < cn(|)^ ^ 
for all j. Finally, we obtain for the total runtime r: 


1=1 1=1 
OO 


acnl - 


h-j-l 


= acn 


h-2 , 

Y.{\ 


= 0{n) as n —>■ cx) 


. 3 1 

acn- 




21 -| 2^2 1-20 


The above does not apply, if 0 does not comply with 
any period length smaller than n — 1. This is because 
assumption and approximation do not hold any¬ 
more as T is always in the right branch and thus it can¬ 
not be considered random which interval is investigated. 
In this case, even if 0 = 0, we have to check the in- 

and thus the 


tervals (f, f ], (f, f ],..., ( \ 

total runtime is: 


r’^max"! 1 

r= ^ (2i-h l)c = -d^c 

2=2 

— ^max ^ 

In particular, we have r = 0{n^) as n —>■ oo, if Tmax oc n. 
A similar approximation can be made, if the shortest 
period length that 0 complies with is close to n. 
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